#Hawkish Biases Replication 2
#Last modified: September 30, 2021
#Josh Kertzer

## This R file contains the code necessary to replicate the analysis from the supplementary appendix. 
## Run Hawkish Biases Replication 1.R first to generate dataframes, etc.
# All of the following analyses were carried out using R version 3.5.2. GUI 1.75 (El Capitan build 7612) on an M1 Pro Macbook Pro running MacOS Monterey; they were also replicated on a 2017 3 Ghz 10-Core Intel Xeon W iMac Pro. 


### Table 2.2: No evidence of treatment spillover effects

trts <- c("lossTrt_pt", "fatalTrt_ib", "chinaTrt_rd")

#Individual analysis
dvs <- c("riskyDV_pt", "intentDV_ib", "supportDV_rd")

b.mat.i <- matrix(NA, 3,3)
s.mat.i <- matrix(NA, 3,3)
for (i in 1:3){
	x <- dat.i[,trts[i]]
	for (j in i:3){
		y <- dat.i[,dvs[j]]
		temp <- lm(y ~ x)
		b.mat.i[i,j] <- coef(temp)[2]
		s.mat.i[i,j] <- sqrt(diag(vcov(temp)))[2]		
	}
}

rownames(b.mat.i) <- c("PT", "IB", "RD")

#Horizontal analysis
dfs <- c("dat.f.1.pt", "dat.f.1.ib", "dat.f.1.rd")
dvs <- c("medDV_pt", "medDV_ib", "medDV_rd")

b.mat.f <- matrix(NA,3,3)
s.mat.f <- matrix(NA,3,3)
for (i in 1:3){
	x <- get(dfs[i])[,trts[i]]
	for (j in i:3){
		temp.df <- get(dfs[j])
		y <- temp.df[,dvs[j]]
		temp <- lm(y ~ x)
		b.mat.f[i,j] <- coef(temp)[2]
		s.mat.f[i,j] <- sqrt(diag(vcov(temp)))[2]		
	}
}

rownames(b.mat.f) <- c("PT", "IB", "RD")

#Hierarchical analysis
dfs <- c("dat.h.1.pt", "dat.h.1.ib", "dat.h.1.rd")
dvs <- c("riskyDV_pt", "intentDV_ib", "supportDV_rd")

b.mat.h <- matrix(NA,3,3)
s.mat.h <- matrix(NA,3,3)
for (i in 1:3){
	x <- get(dfs[i])[,trts[i]]
	for (j in i:3){
		temp.df <- get(dfs[j])
		y <- temp.df[,dvs[j]]
		temp <- lm(y ~ x)
		b.mat.h[i,j] <- coef(temp)[2]
		s.mat.h[i,j] <- sqrt(diag(vcov(temp)))[2]		
	}
}

rownames(b.mat.h) <- c("PT", "IB", "RD")

### Combine outputs for plot (not the most elegant!)
out.b <- rbind(b.mat.i, b.mat.f, b.mat.h)
out.s <- rbind(s.mat.i, s.mat.f, s.mat.h)

stargazer(round(out.b, digits=3))
which(abs(out.b/out.s)>1.96, arr.ind=TRUE) #Use to denote significant effects

##### Table 2.3: Prospect theory experiment

# Individual vs horizontal
dat.i.f.pt <- bind_rows(dat.i, dat.f.1.pt)
dat.i.f.pt$med_DV <- ifelse(is.na(dat.i.f.pt$medDV_pt), dat.i.f.pt$riskyDV_pt, dat.i.f.pt$medDV_pt) #Create common DV for median voter in horizontal, and individuals
dat.i.f.pt$un_DV <- ifelse(is.na(dat.i.f.pt$un_grp_pt), dat.i.f.pt$riskyDV_pt, ifelse(dat.i.f.pt$un_grp_pt==1,dat.i.f.pt$riskyDV_pt, NA))  #Create common DV for unanimous group in horizontal, and individuals
dat.i.f.pt$maj_DV <- ifelse(is.na(dat.i.f.pt$maj_grp_pt), dat.i.f.pt$riskyDV_pt, ifelse(dat.i.f.pt$maj_grp_pt==1,dat.i.f.pt$riskyDV_pt, NA))  #Create common DV for majority group in horizontal, and individuals
dat.i.f.pt$groupCond <- as.factor(dat.i.f.pt$groupCond)
dat.i.f.pt$groupCond <- relevel(dat.i.f.pt$groupCond, ref="Individual")

#Individual vs horizontal (median voter)
mod.if0 <- lm(med_DV ~ lossTrt_pt * groupCond, dat=dat.i.f.pt)
mod.if1 <- lm(med_DV ~ lossTrt_pt * groupCond + age + male + educ1 + income1 + christian + white + pid1 + interest1 + cognition1 + sdo1 + rwa1 + aggression1 + risk1 + mi1 + iso1 + persExtra1 + persAgree1 + persConsc1 + persNeuro1 + persOpen1, dat=dat.i.f.pt)
#Individual vs horizontal (unanimous groups)
mod.if2 <- lm(un_DV ~ lossTrt_pt * groupCond + age + male + educ1 + income1 + christian + white + pid1 + interest1 + cognition1 + sdo1 + rwa1 + aggression1 + risk1 + mi1 + iso1 + persExtra1 + persAgree1 + persConsc1 + persNeuro1 + persOpen1, dat=dat.i.f.pt)
#Individual vs horizontal (majority-rule groups)
mod.if3 <- lm(maj_DV ~ lossTrt_pt * groupCond + age + male + educ1 + income1 + christian + white + pid1 + interest1 + cognition1 + sdo1 + rwa1 + aggression1 + risk1 + mi1 + iso1 + persExtra1 + persAgree1 + persConsc1 + persNeuro1 + persOpen1, dat=dat.i.f.pt)

#Individual vs hierarchical: (leader-level comparison)
dat.i.h.pt <- bind_rows(dat.i, dat.h.1.pt)
dat.i.h.pt$groupCond <- as.factor(dat.i.h.pt$groupCond)
dat.i.h.pt$groupCond <- relevel(dat.i.h.pt$groupCond, ref="Individual")

mod.ih0 <- lm(riskyDV_pt ~ lossTrt_pt * groupCond, dat=dat.i.h.pt)
mod.ih1 <- lm(riskyDV_pt ~ lossTrt_pt * groupCond + age + male + educ1 + income1 + christian + white + pid1 + interest1 + cognition1 + sdo1 + rwa1 + aggression1 + risk1 + mi1 + iso1 + persExtra1 + persAgree1 + persConsc1 + persNeuro1 + persOpen1, dat=dat.i.h.pt)

#Individual vs hierarchical: (group-level comparison)
dat.h.1.pt2 <- dat.h.1.pt
a <- grep("mi1", colnames(dat.h.1.pt2))[1] #First remove leader-level covariates so we can merge
b <- grep("smallgrp1", colnames(dat.h.1.pt2))
dat.h.1.pt2 <- dat.h.1.pt2[,-(a:b)]
loc <- grep("_g", colnames(dat.h.1.pt2))
colnames(dat.h.1.pt2)[loc] <- gsub("_g","",colnames(dat.h.1.pt2)[loc])
dat.i.h.pt2 <- bind_rows(dat.i, dat.h.1.pt2)
dat.i.h.pt2$groupCond <- as.factor(dat.i.h.pt2$groupCond)
dat.i.h.pt2$groupCond <- relevel(dat.i.h.pt2$groupCond, ref="Individual")
mod.ih2 <- lm(riskyDV_pt ~ lossTrt_pt * groupCond + age + male + educ1 + income1 + christian + white + pid1 + interest1 + cognition1 + sdo1 + rwa1 + aggression1 + risk1 + mi1 + iso1 + persExtra1 + persAgree1 + persConsc1 + persNeuro1 + persOpen1, dat=dat.i.h.pt2)

stargazer(mod.if0, mod.if1, mod.if2, mod.if3, mod.ih0, mod.ih1, mod.ih2, title="Prospect theory experiment", style="apsr", digits=3, label="tab:pt", covariate.labels=c("Loss frame", "Horizontal condition", "Age", "Male", "Education", "Income", "Christian", "White", "Party ID", "Political interest", "Cognition", "SDO", "RWA", "Aggression", "Risk attitudes", "MI", "Isolationism", "Extraversion", "Agreeableness", "Conscientiousness", "Neuroticism", "Openness", "Loss x Horizontal", "Hierarchical", "Loss x Hierarchical", "Constant"), dep.var.labels.include=FALSE, column.labels=c("Horizontal", "Hierarchical"), column.separate=c(4,3)) #Hide controls in published version to make the table easier to read

#### Table 2.4: Intentionality bias experiment

# Individual vs horizontal 
dat.i.f.ib <- bind_rows(dat.i, dat.f.1.ib)
dat.i.f.ib$med_DV <- ifelse(is.na(dat.i.f.ib$medDV_ib), dat.i.f.ib$intentDV_ib, dat.i.f.ib$medDV_ib) #Create common DV for median voter in horizontal, and individuals
dat.i.f.ib$un_DV <- ifelse(is.na(dat.i.f.ib$un_grp_ib), dat.i.f.ib$intentDV_ib, ifelse(dat.i.f.ib$un_grp_ib==1,dat.i.f.ib$intentDV_ib, NA))  #Create common DV for unanimous group in horizontal, and individuals
dat.i.f.ib$maj_DV <- ifelse(is.na(dat.i.f.ib$maj_grp_ib), dat.i.f.ib$intentDV_ib, ifelse(dat.i.f.ib$maj_grp_ib==1,dat.i.f.ib$intentDV_ib, NA))  #Create common DV for majority group in horizontal, and individuals
dat.i.f.ib$groupCond <- as.factor(dat.i.f.ib$groupCond)
dat.i.f.ib$groupCond <- relevel(dat.i.f.ib$groupCond, ref="Individual")
#Individual vs horizontal (median voter)
mod.if0 <- lm(med_DV ~ fatalTrt_ib * groupCond, dat=dat.i.f.ib)
mod.if1 <- lm(med_DV ~ fatalTrt_ib * groupCond + age + male + educ1 + income1 + christian + white + pid1 + interest1 + cognition1 + sdo1 + rwa1 + aggression1 + risk1 + mi1 + iso1 + persExtra1 + persAgree1 + persConsc1 + persNeuro1 + persOpen1, dat=dat.i.f.ib)
#Individual vs horizontal (unanimous)
mod.if2 <- lm(un_DV ~ fatalTrt_ib * groupCond + age + male + educ1 + income1 + christian + white + pid1 + interest1 + cognition1 + sdo1 + rwa1 + aggression1 + risk1 + mi1 + iso1 + persExtra1 + persAgree1 + persConsc1 + persNeuro1 + persOpen1, dat=dat.i.f.ib)
#Individual vs horizontal (majority-rule)
mod.if3 <- lm(med_DV ~ fatalTrt_ib * groupCond + age + male + educ1 + income1 + christian + white + pid1 + interest1 + cognition1 + sdo1 + rwa1 + aggression1 + risk1 + mi1 + iso1 + persExtra1 + persAgree1 + persConsc1 + persNeuro1 + persOpen1, dat=dat.i.f.ib)

#Individual vs hierarchical: (leader-level comparison)
dat.i.h.ib <- bind_rows(dat.i, dat.h.1.ib)
dat.i.h.ib$groupCond <- as.factor(dat.i.h.ib$groupCond)
dat.i.h.ib$groupCond <- relevel(dat.i.h.ib$groupCond, ref="Individual")

mod.ih0 <- lm(intentDV_ib ~ fatalTrt_ib * groupCond, dat=dat.i.h.ib)
mod.ih1 <- lm(intentDV_ib ~ fatalTrt_ib * groupCond + age + male + educ1 + income1 + christian + white + pid1 + interest1 + cognition1 + sdo1 + rwa1 + aggression1 + risk1 + mi1 + iso1 + persExtra1 + persAgree1 + persConsc1 + persNeuro1 + persOpen1, dat=dat.i.h.ib)

#Individual vs hierarchical: (group-level comparison)
dat.h.1.ib2 <- dat.h.1.ib
a <- grep("mi1", colnames(dat.h.1.ib2))[1] #First remove leader-level covariates so we can merge
b <- grep("smallgrp1", colnames(dat.h.1.ib2))
dat.h.1.ib2 <- dat.h.1.ib2[,-(a:b)]
loc <- grep("_g", colnames(dat.h.1.ib2))
colnames(dat.h.1.ib2)[loc] <- gsub("_g","",colnames(dat.h.1.ib2)[loc])
dat.i.h.ib2 <- bind_rows(dat.i, dat.h.1.ib2)
dat.i.h.ib2$groupCond <- as.factor(dat.i.h.ib2$groupCond)
dat.i.h.ib2$groupCond <- relevel(dat.i.h.ib2$groupCond, ref="Individual")

mod.ih2 <- lm(intentDV_ib ~ fatalTrt_ib * groupCond + age + male + educ1 + income1 + christian + white + pid1 + interest1 + cognition1 + sdo1 + rwa1 + aggression1 + risk1 + mi1 + iso1 + persExtra1 + persAgree1 + persConsc1 + persNeuro1 + persOpen1, dat=dat.i.h.ib2)

stargazer(mod.if0, mod.if1, mod.if2, mod.if3, mod.ih0, mod.ih1, mod.ih2, title="Intentionality bias experiment", style="apsr", digits=3, label="tab:ib", covariate.labels=c("Fatalities", "Horizontal condition", "Age", "Male", "Education", "Income", "Christian", "White", "Party ID", "Political interest", "Cognition", "SDO", "RWA", "Aggression", "Risk attitudes", "MI", "Isolationism", "Extraversion", "Agreeableness", "Conscientiousness", "Neuroticism", "Openness", "Fatalities x Horizontal", "Hierarchical", "Fatalities x Hierarchical", "Constant"), dep.var.labels.include=FALSE, column.labels=c("Horizontal", "Hierarchical"), column.separate=c(4,3)) #Hide controls in published version to make the table easier to read

#### Table 2.5: Reactive devaluation experiment

# Individual vs horizontal
dat.i.f.rd <- bind_rows(dat.i, dat.f.1.rd)
dat.i.f.rd$med_DV <- ifelse(is.na(dat.i.f.rd$medDV_rd), dat.i.f.rd$supportDV_rd, dat.i.f.rd$medDV_rd) #Create common DV for median voter in horizontal, and individuals
dat.i.f.rd$un_DV <- ifelse(is.na(dat.i.f.rd$un_grp_rd), dat.i.f.rd$supportDV_rd, ifelse(dat.i.f.rd$un_grp_rd==1,dat.i.f.rd$supportDV_rd, NA))  #Create common DV for unanimous group in horizontal, and individuals
dat.i.f.rd$maj_DV <- ifelse(is.na(dat.i.f.rd$maj_grp_rd), dat.i.f.rd$supportDV_rd, ifelse(dat.i.f.rd$maj_grp_rd==1,dat.i.f.rd$supportDV_rd, NA))  #Create common DV for majority group in horizontal, and individuals
dat.i.f.rd$groupCond <- as.factor(dat.i.f.rd$groupCond)
dat.i.f.rd$groupCond <- relevel(dat.i.f.rd$groupCond, ref="Individual")
#Individual vs horizontal (median voter)
mod.if0 <- lm(med_DV ~ chinaTrt_rd * groupCond, dat=dat.i.f.rd)
mod.if1 <- lm(med_DV ~ chinaTrt_rd * groupCond + age + male + educ1 + income1 + christian + white + pid1 + interest1 + cognition1 + sdo1 + rwa1 + aggression1 + risk1 + mi1 + iso1 + persExtra1 + persAgree1 + persConsc1 + persNeuro1 + persOpen1, dat=dat.i.f.rd)
#Individual vs horizontal (unanimous)
mod.if2 <- lm(un_DV ~ chinaTrt_rd * groupCond + age + male + educ1 + income1 + christian + white + pid1 + interest1 + cognition1 + sdo1 + rwa1 + aggression1 + risk1 + mi1 + iso1 + persExtra1 + persAgree1 + persConsc1 + persNeuro1 + persOpen1, dat=dat.i.f.rd)
#Individual vs horizontal (majority-rule)
mod.if3 <- lm(maj_DV ~ chinaTrt_rd * groupCond + age + male + educ1 + income1 + christian + white + pid1 + interest1 + cognition1 + sdo1 + rwa1 + aggression1 + risk1 + mi1 + iso1 + persExtra1 + persAgree1 + persConsc1 + persNeuro1 + persOpen1, dat=dat.i.f.rd)
#Individual vs hierarchical: (leader-level comparison)
dat.i.h.rd <- bind_rows(dat.i, dat.h.1.rd)
dat.i.h.rd$groupCond <- as.factor(dat.i.h.rd$groupCond)
dat.i.h.rd$groupCond <- relevel(dat.i.h.rd$groupCond, ref="Individual")
mod.ih0 <- lm(supportDV_rd ~ chinaTrt_rd * groupCond, dat=dat.i.h.rd)
mod.ih1 <- lm(supportDV_rd ~ chinaTrt_rd * groupCond + age + male + educ1 + income1 + christian + white + pid1 + interest1 + cognition1 + sdo1 + rwa1 + aggression1 + risk1 + mi1 + iso1 + persExtra1 + persAgree1 + persConsc1 + persNeuro1 + persOpen1, dat=dat.i.h.rd)

#Individual vs hierarchical: (group-level comparison)
dat.h.1.rd2 <- dat.h.1.rd
a <- grep("mi1", colnames(dat.h.1.rd2))[1] #First remove leader-level covariates so we can merge
b <- grep("smallgrp1", colnames(dat.h.1.rd2))
dat.h.1.rd2 <- dat.h.1.rd2[,-(a:b)]
loc <- grep("", colnames(dat.h.1.rd2))
colnames(dat.h.1.rd2)[loc] <- gsub("_g","",colnames(dat.h.1.rd2)[loc])
dat.i.h.rd2 <- bind_rows(dat.i, dat.h.1.rd2)
dat.i.h.rd2 <- bind_rows(dat.i, dat.h.1.rd2)
dat.i.h.rd2$groupCond <- as.factor(dat.i.h.rd2$groupCond)
dat.i.h.rd2$groupCond <- relevel(dat.i.h.rd2$groupCond, ref="Individual")
mod.ih2 <- lm(supportDV_rd ~ chinaTrt_rd * groupCond + age + male + educ1 + income1 + christian + white + pid1 + interest1 + cognition1 + sdo1 + rwa1 + aggression1 + risk1 + mi1 + iso1 + persExtra1 + persAgree1 + persConsc1 + persNeuro1 + persOpen1, dat=dat.i.h.rd2)

stargazer(mod.if0, mod.if1, mod.if2, mod.if3, mod.ih0, mod.ih1, mod.ih2, title="Reactive devaluation experiment", style="apsr", digits=3, label="tab:ib", covariate.labels=c("China authored", "Horizontal condition", "Age", "Male", "Education", "Income", "Christian", "White", "Party ID", "Political interest", "Cognition", "SDO", "RWA", "Aggression", "Risk attitudes", "MI", "Isolationism", "Extraversion", "Agreeableness", "Conscientiousness", "Neuroticism", "Openness", "China x Horizontal", "Hierarchical", "China x Hierarchical", "Constant"), dep.var.labels.include=FALSE, column.labels=c("Horizontal", "Hierarchical"), column.separate=c(4,3)) #Hide controls in published version to make the table easier to read

###### Sensitivity analysis
#Load and generate alternate datasets
dat.c00 <- get(load("Hawkish-Biases-2.RData"))

dat.i.00 <- dat.c00[which(dat.c00$groupCond=="Individual"),] #Individual condition
dat.f.00 <- dat.c00[which(dat.c00$groupCond=="Horizontal"),] #Horizontal condition
dat.h.00 <- dat.c00[which(dat.c00$groupCond=="Hierarchical"),] #Hierarchical condition

#Hierarchical
dat.h.00 <- dat.c00[which(dat.c00$groupCond=="Hierarchical" & dat.c00$role_pt=="Leader"),]

#Horizontal
dat.f.00.pt <- data.frame(GroupID_pt=unique(dat.f.00$GroupID_pt), groupCond="Horizontal")
dat.f.00.pt$medDV_pt <- NA #Median voter
dat.f.00.pt$lossTrt_pt <- NA #Treatment

for (i in unique(dat.f.00$GroupID_pt)){
	j <- which(dat.f.00$GroupID_pt==i)
	k <- which(dat.f.00.pt$GroupID_pt==i)
	group.size <- length(na.omit(dat.f.00$riskyDV_pt[j]))
	dat.f.00.pt$groupN_pt[k] <- group.size	
	dat.f.00.pt$medDV_pt[k] <- median(dat.f.00$riskyDV_pt[j], na.rm=TRUE)
	dat.f.00.pt$gainTrt_pt[k] <- ifelse(length(unique(dat.f.00$gainTrt_pt[j]))==1,unique(dat.f.00$gainTrt_pt[j]),NA)
	dat.f.00.pt$lossTrt_pt[k] <- ifelse(length(unique(dat.f.00$lossTrt_pt[j]))==1,unique(dat.f.00$lossTrt_pt[j]),NA)		
}

dat.f.00.ib <- data.frame(GroupID_ib=unique(dat.f.00$GroupID_ib), groupCond="Horizontal")
dat.f.00.ib$medDV_ib <- NA #Median voter
dat.f.00.ib$fatalTrt_ib <- NA #Treatments
dat.f.00.ib$groupN_ib <- NA #Group size

for (i in unique(dat.f.00$GroupID_ib)){
	j <- which(dat.f.00$GroupID_ib==i)
	k <- which(dat.f.00.ib$GroupID_ib==i)
	group.size <- length(na.omit(dat.f.00$intentDV_ib[j]))
	dat.f.00.ib$groupN_ib[k] <- group.size
	dat.f.00.ib$medDV_ib[k] <- median(dat.f.00$intentDV_ib[j], na.rm=TRUE)		
	dat.f.00.ib$fatalTrt_ib[k] <- ifelse(length(unique(dat.f.00$fatalTrt_ib[j]))==1,unique(dat.f.00$fatalTrt_ib[j]),NA)				
}

#RD submodule	
dat.f.00.rd <- data.frame(GroupID_rd=unique(dat.f.00$GroupID_rd), groupCond="Horizontal")
dat.f.00.rd$medDV_rd <- NA #Median voter
dat.f.00.rd$chinaTrt_rd <- NA #Treatments
dat.f.00.rd$groupN_rd <- NA #Group size

for (i in unique(dat.f.00$GroupID_rd)){
	j <- which(dat.f.00$GroupID_rd==i)
	k <- which(dat.f.00.rd$GroupID_rd==i)
	group.size <- length(na.omit(dat.f.00$supportDV_rd[j]))
	dat.f.00.rd$groupN_rd[k] <- group.size
	dat.f.00.rd$medDV_rd[k] <- median(dat.f.00$supportDV_rd[j],na.rm=TRUE)			
	dat.f.00.rd$chinaTrt_rd[k] <- ifelse(length(unique(dat.f.00$chinaTrt_rd[j]))==1,unique(dat.f.00$chinaTrt_rd[j]),NA)
}

#Now generate extreme bound equivalents

#PT submodule
dat.f.00.pt$medDV0_pt <- NA #Median vote (where missing values imputed as 0)
dat.f.00.pt$medDV1_pt <- NA #Median vote (where missing values imputed as 1)
dat.i$riskyDV0_pt <- ifelse(is.na(dat.i$riskyDV_pt),0,dat.i$riskyDV_pt) #missing values imputed as 0
dat.i$riskyDV1_pt <- ifelse(is.na(dat.i$riskyDV_pt),1,dat.i$riskyDV_pt) #missing values imputed as 1
dat.h.00$riskyDV0_pt <- ifelse(is.na(dat.h.00$riskyDV_pt),0,dat.h.00$riskyDV_pt) #missing leader values imputed as 0
dat.h.00$riskyDV1_pt <- ifelse(is.na(dat.h.00$riskyDV_pt),1,dat.h.00$riskyDV_pt) #missing leader values imputed as 1

for (i in unique(dat.f.00$GroupID_pt)){
	j <- which(dat.f.00$GroupID_pt==i)
	k <- which(dat.f.00.pt$GroupID_pt==i)
	temp <- dat.f.00$riskyDV_pt[j]
	dat.f.00.pt$medDV0_pt[k] <- median(ifelse(is.na(temp),0,temp))
	dat.f.00.pt$medDV1_pt[k] <- median(ifelse(is.na(temp),1,temp))
}

#IB submodule
dat.f.00.ib$medDV0_ib <- NA #Median vote (where missing values imputed as 0)
dat.f.00.ib$medDV1_ib <- NA #Median vote (where missing values imputed as 1)
dat.i$intentDV0_ib <- ifelse(is.na(dat.i$intentDV_ib),0,dat.i$intentDV_ib) #missing values imputed as 0
dat.i$intentDV1_ib <- ifelse(is.na(dat.i$intentDV_ib),1,dat.i$intentDV_ib) #missing values imputed as 1
dat.h.00$intentDV0_ib <- ifelse(is.na(dat.h.00$intentDV_ib),0,dat.h.00$intentDV_ib) #missing leader values imputed as 0
dat.h.00$intentDV1_ib <- ifelse(is.na(dat.h.00$intentDV_ib),1,dat.h.00$intentDV_ib) #missing leader values imputed as 1

for (i in unique(dat.f.00$GroupID_ib)){
	j <- which(dat.f.00$GroupID_ib==i)
	k <- which(dat.f.00.ib$GroupID_ib==i)
	temp <- dat.f.00$intentDV_ib[j]
	dat.f.00.ib$medDV0_ib[k] <- median(ifelse(is.na(temp),0,temp))
	dat.f.00.ib$medDV1_ib[k] <- median(ifelse(is.na(temp),1,temp))
}

#RD submodule	
dat.f.00.rd$medDV0_rd <- NA #Median vote (where missing values imputed as 0)
dat.f.00.rd$medDV1_rd <- NA #Median vote (where missing values imputed as 1)
dat.i$supportDV0_rd <- ifelse(is.na(dat.i$supportDV_rd),0,dat.i$supportDV_rd) #missing values imputed as 0
dat.i$supportDV1_rd <- ifelse(is.na(dat.i$supportDV_rd),1,dat.i$supportDV_rd) #missing values imputed as 1
dat.h.00$supportDV0_rd <- ifelse(is.na(dat.h.00$supportDV_rd),0,dat.h.00$supportDV_rd) #missing leader values imputed as 0
dat.h.00$supportDV1_rd <- ifelse(is.na(dat.h.00$supportDV_rd),1,dat.h.00$supportDV_rd) #missing leader values imputed as 1

for (i in unique(dat.f.00$GroupID_rd)){
	j <- which(dat.f.00$GroupID_rd==i)
	k <- which(dat.f.00.rd$GroupID_rd==i)
	temp <- dat.f.00$supportDV_rd[j]
	dat.f.00.rd$medDV0_rd[k] <- median(ifelse(is.na(temp),0,temp))
	dat.f.00.rd$medDV1_rd[k] <- median(ifelse(is.na(temp),1,temp))
}

#Simulate counterfactual: given what we know about the relationship between respondent demographics and their choice in the individual condition, what would these respondents do in the group conditions assuming they behaved similarly?


#Generate slim version of the data with no rows of complete missingness
ind.pt <- dat.i %>% dplyr::select(riskyDV_pt, lossTrt_pt, mi1, iso1, risk1, sdo1, rwa1, persExtra1, persAgree1, persConsc1, persNeuro1, persOpen1, cognition1, aggression1, male, age, educ1, income1, white, pid1, ideo1, interest1) %>% drop_na(everything())
ind.pt$riskyDV_pt <- as.factor(ind.pt$riskyDV_pt) #Generate a factor version for classification

ind.ib <- dat.i %>% dplyr::select(intentDV_ib, fatalTrt_ib, mi1, iso1, risk1, sdo1, rwa1, persExtra1, persAgree1, persConsc1, persNeuro1, persOpen1, cognition1, aggression1, male, age, educ1, income1, white, pid1, ideo1, interest1) %>% drop_na(everything())

ind.rd <- dat.i %>% dplyr::select(supportDV_rd, chinaTrt_rd, mi1, iso1, risk1, sdo1, rwa1, persExtra1, persAgree1, persConsc1, persNeuro1, persOpen1, cognition1, aggression1, male, age, educ1, income1, white, pid1, ideo1, interest1) %>% drop_na(everything())

## For PT, use a neural net (but for classification)
set.seed(43215)
test.grid.nn <- expand.grid(size=c(1,2,3,4,5,6,7,8,9,10,11,12), decay=c(2,1,0.1,0.01,0.001))
ind.pt.nn <- train(riskyDV_pt ~ ., data=ind.pt, method="nnet", na.action=na.omit, tuneGrid=test.grid.nn, trControl=trainControl(method="repeatedcv", repeats=5))
ind.pt.nn #size=10, decay = 1, using accuracy criterion (0.66)

nn.f.pred.pt <- predict(ind.pt.nn, newdata=dat.f.00) #Predicted DV for horizontal condition
nn.f.pred.pt <- as.numeric(nn.f.pred.pt)-1 #Convert from factor 
cor(nn.f.pred.pt, dat.f.00$riskyDV_pt, use="complete.obs", method="spearman") #rho=0.29

nn.h.pred.pt <- predict(ind.pt.nn, newdata=dat.h.00) #Predicted DV for hierarchical condition
nn.h.pred.pt <- as.numeric(nn.h.pred.pt)-1 #Convert from factor 
cor(nn.h.pred.pt, dat.h.00$riskyDV_pt, use="complete.obs", method="spearman") #rho=0.37

## For IB, use a neural net
set.seed(43215)
ind.nn <- train(intentDV_ib ~ ., data=ind.ib, method="nnet", na.action=na.omit, tuneGrid=test.grid.nn, trControl=trainControl(method="repeatedcv", repeats=5))
ind.nn #size=4, decay = 0.1

nn.f.pred.ib <- predict(ind.nn, newdata=dat.f.00) #Predicted DV for horizontal condition
cor(nn.f.pred.ib, dat.f.00$intentDV_ib, use="complete.obs") #r=0.22
nn.h.pred.ib <- predict(ind.nn, newdata=dat.h.00) #Predicted DV for hierarchical condition
cor(nn.h.pred.ib, dat.h.00$intentDV_ib, use="complete.obs") #r=0.22

## For RD, use a neural net
set.seed(43215)
ind.rd.nn <- train(supportDV_rd ~ ., data=ind.rd, method="nnet", na.action=na.omit, tuneGrid=test.grid.nn, trControl=trainControl(method="repeatedcv", repeats=5))
ind.rd.nn #size=6, decay = 0.1

nn.f.pred.rd <- predict(ind.rd.nn, newdata=dat.f.00) #Predicted DV for horizontal condition
cor(nn.f.pred.rd, dat.f.00$supportDV_rd, use="complete.obs") #r=0.18
nn.h.pred.rd <- predict(ind.rd.nn, newdata=dat.h.00) #Predicted DV for horizontal condition
cor(nn.h.pred.rd, dat.h.00$supportDV_rd, use="complete.obs") #r=0.13

## Now, merge these into missing observations
i <- which(is.na(dat.f.00$riskyDV_pt))
dat.f.00$riskyDVnn_pt <- NA
dat.f.00$riskyDVnn_pt[i] <- nn.f.pred.pt[i]

i <- which(is.na(dat.h.00$riskyDV_pt))
dat.h.00$riskyDVnn_pt <- NA
dat.h.00$riskyDVnn_pt[i] <- nn.h.pred.pt[i]

i <- which(is.na(dat.f.00$intentDV_ib))
dat.f.00$intentDVnn_ib <- NA
dat.f.00$intentDVnn_ib[i] <- nn.f.pred.ib[i]

i <- which(is.na(dat.h.00$intentDV_ib))
dat.h.00$intentDVnn_ib <- NA
dat.h.00$intentDVnn_ib[i] <- nn.h.pred.ib[i]

i <- which(is.na(dat.f.00$supportDV_rd))
dat.f.00$supportDVnn_rd <- NA
dat.f.00$supportDVnn_rd[i] <- nn.f.pred.rd[i]

i <- which(is.na(dat.h.00$supportDV_rd))
dat.h.00$supportDVnn_rd <- NA
dat.h.00$supportDVnn_rd[i] <- nn.h.pred.rd[i]

## Now, for horizontal, two different scenarios: a) antisocial case: these respondents behave as if they were individuals in the group conditions; b) influencer case: these respondents were the pivotal member of the group (only applies in horizontal; this is automatically the case for hierarchical)

#PT submodule
dat.f.00.pt$medDVA_pt <- NA #Median vote (scenario A)
dat.f.00.pt$medDVB_pt <- NA #Median vote (scenario B)
dat.h.00$riskyDVAB_pt <- ifelse(is.na(dat.h.00$riskyDV_pt),dat.h.00$riskyDVnn_pt,dat.h.00$riskyDV_pt) #Leaders behave as they would in individual

for (i in unique(dat.f.00$GroupID_pt)){
	j <- which(dat.f.00$GroupID_pt==i)
	k <- which(dat.f.00.pt$GroupID_pt==i)
	temp <- dat.f.00$riskyDV_pt[j]
	if(any(is.na(temp))){
		dat.f.00.pt$medDVA_pt[k] <- median(c(temp,dat.f.00$riskyDVnn_pt[j]),na.rm=TRUE) #Incorporate NN prediction into median vote
		dat.f.00.pt$medDVB_pt[k] <- median(dat.f.00$riskyDVnn_pt[j],na.rm=TRUE) #Make NN prediction group decision (take median of NN predictions, in case of multiple missing members)
	} else {
		dat.f.00.pt$medDVA_pt[k] <- median(temp)
		dat.f.00.pt$medDVB_pt[k] <- median(temp)
	}
}

#IB submodule
dat.f.00.ib$medDVA_ib <- NA #Median vote (scenario A)
dat.f.00.ib$medDVB_ib <- NA #Median vote (scenario B)
dat.h.00$intentDVAB_ib <- ifelse(is.na(dat.h.00$intentDV_ib),dat.h.00$intentDVnn_ib,dat.h.00$intentDV_ib) #Leaders behave as they would in individual

for (i in unique(dat.f.00$GroupID_ib)){
	j <- which(dat.f.00$GroupID_ib==i)
	k <- which(dat.f.00.ib$GroupID_ib==i)
	temp <- dat.f.00$intentDV_ib[j]
	if(any(is.na(temp))){
		dat.f.00.ib$medDVA_ib[k] <- median(c(temp,dat.f.00$intentDVnn_ib[j]),na.rm=TRUE) #Incorporate NN prediction into median vote
		dat.f.00.ib$medDVB_ib[k] <- median(dat.f.00$intentDVnn_ib[j],na.rm=TRUE) #Make NN prediction group decision (take median of NN predictions, in case of multiple missing members)
	} else {
		dat.f.00.ib$medDVA_ib[k] <- median(temp)
		dat.f.00.ib$medDVB_ib[k] <- median(temp)
	}
}

#RD submodule
dat.f.00.rd$medDVA_rd <- NA #Median vote (scenario A)
dat.f.00.rd$medDVB_rd <- NA #Median vote (scenario B)
dat.h.00$supportDVAB_rd <- ifelse(is.na(dat.h.00$supportDV_rd),dat.h.00$supportDVnn_rd,dat.h.00$supportDV_rd) #Leaders behave as they would in individual

for (i in unique(dat.f.00$GroupID_rd)){
	j <- which(dat.f.00$GroupID_rd==i)
	k <- which(dat.f.00.rd$GroupID_rd==i)
	temp <- dat.f.00$supportDV_rd[j]
	if(any(is.na(temp))){
		dat.f.00.rd$medDVA_rd[k] <- median(c(temp,dat.f.00$supportDVnn_rd[j]),na.rm=TRUE) #Incorporate NN prediction into median vote
		dat.f.00.rd$medDVB_rd[k] <- median(dat.f.00$supportDVnn_rd[j],na.rm=TRUE) #Make NN prediction group decision (take median of NN predictions, in case of multiple missing members)
	} else {
		dat.f.00.rd$medDVA_rd[k] <- median(temp)
		dat.f.00.rd$medDVB_rd[k] <- median(temp)
	}
}

### Table 2.6 Prospect theory sensitivity analysis

pt.i <- lm(riskyDV_pt ~ lossTrt_pt, data=dat.i)
pt.im <- lm(riskyDV_pt ~ lossTrt_pt, data=dat.i.00)
pt.i0 <- lm(riskyDV0_pt ~ lossTrt_pt, data=dat.i) 
pt.i1 <- lm(riskyDV1_pt ~ lossTrt_pt, data=dat.i)
pt.f <- lm(medDV_pt ~ lossTrt_pt, data=dat.f.1.pt)
pt.fm <- lm(medDV_pt ~ lossTrt_pt, data=dat.f.00.pt)
pt.f00.0 <- lm(medDV0_pt ~ lossTrt_pt, data=dat.f.00.pt)
pt.f00.1 <- lm(medDV1_pt ~ lossTrt_pt, data=dat.f.00.pt) 
pt.f00.A <- lm(medDVA_pt ~ lossTrt_pt, data=dat.f.00.pt) #X
pt.f00.B <- lm(medDVA_pt ~ lossTrt_pt, data=dat.f.00.pt) #X
pt.h <- lm(riskyDV_pt ~ lossTrt_pt, data=dat.h.1)
pt.hm <- lm(riskyDV_pt ~ lossTrt_pt, data=dat.h.00)
pt.h00.0 <- lm(riskyDV0_pt ~ lossTrt_pt, data=dat.h.00) 
pt.h00.1 <- lm(riskyDV1_pt ~ lossTrt_pt, data=dat.h.00)
pt.h00.AB <- lm(riskyDVAB_pt ~ lossTrt_pt, data=dat.h.00) #X

#Full table
stargazer(pt.i, pt.i0, pt.i1, pt.f, pt.fm, pt.f00.0, pt.f00.1, pt.f00.A, pt.f00.B, pt.h,pt.hm, pt.h00.0, pt.h00.1, pt.h00.AB, title="Prospect theory sensitivity analysis", style="apsr", digits=3, dep.var.labels.include=FALSE, covariate.labels=c("Loss Trt", "Intercept"),column.labels=c("Individual", "Horizontal", "Hierarchical"), omit.stat=c("LL", "ser", "f"),column.separate=c(3,6,5), add.lines=list(c("Bounds"," ","K0","K1","","M","K0","K1","NNA","NNB","","M","K0","K1","NNAB")))

#Reformat to fit on a page
stargazer(pt.i, pt.i0, pt.i1, pt.f, pt.fm, pt.f00.0, pt.f00.1, title="Prospect theory sensitivity analysis", style="apsr", digits=3, dep.var.labels.include=FALSE, covariate.labels=c("Loss Trt", "Intercept"),column.labels=c("Individual", "Horizontal"), omit.stat=c("LL", "ser", "f"),column.separate=c(3,4), add.lines=list(c("Model","Base","Extreme Min","Extreme Max","Base","Complete Obs","Extreme Min","Extreme Max")))
stargazer(pt.f00.A, pt.f00.B, pt.h,pt.hm, pt.h00.0, pt.h00.1, pt.h00.AB, style="apsr", digits=3, dep.var.labels.include=FALSE, covariate.labels=c("Loss Trt", "Intercept"),column.labels=c("Horizontal", "Hierarchical"), omit.stat=c("LL", "ser", "f"),column.separate=c(2,5), add.lines=list(c("Model","Neural Net Antisocial","Neural Net Influencer","Base","Complete Obs","Extreme Min","Extreme Max","Neural net")))



### Table 2.7 Intentionality bias sensitivity analysis

ib.i <- lm(intentDV_ib ~ fatalTrt_ib, data=dat.i)
ib.i0 <- lm(intentDV0_ib ~ fatalTrt_ib, data=dat.i)
ib.i1 <- lm(intentDV1_ib ~ fatalTrt_ib, data=dat.i)
ib.f <- lm(medDV_ib ~ fatalTrt_ib, data=dat.f.1.ib)
ib.fm <- lm(medDV_ib ~ fatalTrt_ib, data=dat.f.00.ib)
ib.f00.0 <- lm(medDV0_ib ~ fatalTrt_ib, data=dat.f.00.ib)
ib.f00.1 <- lm(medDV1_ib ~ fatalTrt_ib, data=dat.f.00.ib)
ib.f00.A <- lm(medDVA_ib ~ fatalTrt_ib, data=dat.f.00.ib)
ib.f00.B <- lm(medDVB_ib ~ fatalTrt_ib, data=dat.f.00.ib)
ib.h <- lm(intentDV_ib ~ fatalTrt_ib, data=dat.h.1)
ib.hm <- lm(intentDV_ib ~ fatalTrt_ib, data=dat.h.00)
ib.h00.0 <- lm(intentDV0_ib ~ fatalTrt_ib, data=dat.h.00)
ib.h00.1 <- lm(intentDV1_ib ~ fatalTrt_ib, data=dat.h.00)
ib.h00.AB <- lm(intentDVAB_ib ~ fatalTrt_ib, data=dat.h.00)


#Full table
stargazer(ib.i, ib.i0, ib.i1, ib.f, ib.fm, ib.f00.0, ib.f00.1, ib.f00.A, ib.f00.B, ib.h,ib.hm, ib.h00.0, ib.h00.1, ib.h00.AB, title="Intentionality bias sensitivity analysis", style="apsr", digits=3, dep.var.labels.include=FALSE, covariate.labels=c("Fatal Trt", "Intercept"),column.labels=c("Individual", "Horizontal", "Hierarchical"), omit.stat=c("LL", "ser", "f"),column.separate=c(3,5,4), add.lines=list(c("Bounds"," ","K0","K1","","K0","K1","NNA","NNB","","K0","K1","NNAB")))
#In pieces to fit on the page
stargazer(ib.i, ib.i0, ib.i1, ib.f, ib.fm, ib.f00.0, ib.f00.1, title="Intentionality bias sensitivity analysis",style="apsr", digits=3, dep.var.labels.include=FALSE, covariate.labels=c("Fatal Trt", "Intercept"),column.labels=c("Individual", "Horizontal"), omit.stat=c("LL", "ser", "f"),column.separate=c(3,4), add.lines=list(c("Model","Base","Extreme Min","Extreme Max","Base","Complete Obs","Extreme Min","Extreme Max")))
stargazer(ib.f00.A, ib.f00.B, ib.h,ib.hm, ib.h00.0, ib.h00.1, ib.h00.AB, style="apsr", digits=3, dep.var.labels.include=FALSE, covariate.labels=c("Fatal Trt", "Intercept"),column.labels=c("Horizontal", "Hierarchical"), omit.stat=c("LL", "ser", "f"),column.separate=c(2,5), add.lines=list(c("Model","Neural Net Antisocial","Neural Net Influencer","Base","Complete Obs","Extreme Min","Extreme Max","Neural net")))


### Table 2.8 Reactive devaluation sensitivity analysis

rd.i <- lm(supportDV_rd ~ chinaTrt_rd, data=dat.i)
rd.i0 <- lm(supportDV0_rd ~ chinaTrt_rd, data=dat.i)
rd.i1 <- lm(supportDV1_rd ~ chinaTrt_rd, data=dat.i)
rd.f <- lm(medDV_rd ~ chinaTrt_rd, data=dat.f.1.rd)
rd.fm <- lm(medDV_rd ~ chinaTrt_rd, data=dat.f.00.rd)
rd.f00.0 <- lm(medDV0_rd ~ chinaTrt_rd, data=dat.f.00.rd)
rd.f00.1 <- lm(medDV1_rd ~ chinaTrt_rd, data=dat.f.00.rd)
rd.f00.A <- lm(medDVA_rd ~ chinaTrt_rd, data=dat.f.00.rd)
rd.f00.B <- lm(medDVB_rd ~ chinaTrt_rd, data=dat.f.00.rd)
rd.h <- lm(supportDV_rd ~ chinaTrt_rd, data=dat.h.1)
rd.hm <- lm(supportDV_rd ~ chinaTrt_rd, data=dat.h.00)
rd.h00.0 <- lm(supportDV0_rd ~ chinaTrt_rd, data=dat.h.00)
rd.h00.1 <- lm(supportDV1_rd ~ chinaTrt_rd, data=dat.h.00)
rd.h00.AB <- lm(supportDVAB_rd ~ chinaTrt_rd, data=dat.h.00)

#Full table
stargazer(rd.i, rd.i0, rd.i1, rd.f, rd.fm, rd.f00.0, rd.f00.1, rd.f00.A, rd.f00.B, rd.h,rd.hm, rd.h00.0, rd.h00.1, rd.h00.AB, title="Reactive devaluation sensitivity analysis", style="apsr", digits=3, dep.var.labels.include=FALSE, covariate.labels=c("China Trt", "Intercept"),column.labels=c("Individual", "Horizontal", "Hierarchical"), omit.stat=c("LL", "ser", "f"),column.separate=c(3,5,4), add.lines=list(c("Bounds"," ","K0","K1","","K0","K1","NNA","NNB","","K0","K1","NNAB")))
#In pieces
stargazer(rd.i, rd.i0, rd.i1, rd.f, rd.fm, rd.f00.0, rd.f00.1, title="Reactive devaluation sensitivity analysis",style="apsr", digits=3, dep.var.labels.include=FALSE, covariate.labels=c("China Trt", "Intercept"),column.labels=c("Individual", "Horizontal"), omit.stat=c("LL", "ser", "f"),column.separate=c(3,4), add.lines=list(c("Model","Base","Extreme Min","Extreme Max","Base","Complete Obs","Extreme Min","Extreme Max")))
stargazer(rd.f00.A, rd.f00.B, rd.h,rd.hm, rd.h00.0, rd.h00.1, rd.h00.AB, style="apsr", digits=3, dep.var.labels.include=FALSE, covariate.labels=c("China Trt", "Intercept"),column.labels=c("Horizontal", "Hierarchical"), omit.stat=c("LL", "ser", "f"),column.separate=c(2,5), add.lines=list(c("Model","Neural Net Antisocial","Neural Net Influencer","Base","Complete Obs","Extreme Min","Extreme Max","Neural net")))


### Figure 3.2: Hawkish biases do not aggregate away in larger groups

#Prospect theory group size analysis for individual (generate 1000 groups with replacement (all group members must be in the same treatment condition), each of size N (where N=1-20); calcualte ATE, and resample 1000 times; repeat for each group size)
set.seed(43215)
B <- 1500
ig <- matrix(NA, nrow=B, ncol=20)
ctrl.df <- dat.i[which(dat.i$lossTrt_pt==0),] #Subset of observations in control condition
trt.df <- dat.i[which(dat.i$lossTrt_pt==1),] #Subset of observations in treatment condition
for (h in 1:B){
	for (i in 1:20){
		temp.df <- data.frame(trt=rbinom(1000,1,0.5), DV=NA)	
		for (j in 1:length(temp.df$trt)){
			temp.df$DV[j] <- ifelse(temp.df$trt[j]==0,median(ctrl.df$riskyDV_pt[sample(1:nrow(ctrl.df),i,replace=TRUE)],na.rm=TRUE), median(trt.df$riskyDV_pt[sample(1:nrow(trt.df),i,replace=TRUE)], na.rm=TRUE))
		}
		ig[h,i] <- mean(temp.df$DV[which(temp.df$trt==1)],na.rm=TRUE) - mean(temp.df$DV[which(temp.df$trt==0)],na.rm=TRUE)
	}
}
	
ig2 <- as.data.frame(ig)
colnames(ig2) <- c(1:20)
ig2$id <- 1:nrow(ig2)

ig.df <- pivot_longer(ig2, cols=1:20)
ig.df$name <- as.numeric(ig.df$name)


#Intentionality bias group size analysis for individual (generate 1000 groups with replacement (all group members must be in the same treatment condition), each of size N (where N=1-20); calcualte ATE, and resample 1000 times; repeat for each group size)
set.seed(43215)
B <- 1500
ig.ib <- matrix(NA, nrow=B, ncol=20)
ctrl.df <- dat.i[which(dat.i$fatalTrt_ib==0),] #Subset of observations in control condition
trt.df <- dat.i[which(dat.i$fatalTrt_ib==1),] #Subset of observations in treatment condition
for (h in 1:B){
	for (i in 1:20){
		temp.df <- data.frame(trt=rbinom(1000,1,0.5), DV=NA)	
		for (j in 1:length(temp.df$trt)){
			temp.df$DV[j] <- ifelse(temp.df$trt[j]==0,median(ctrl.df$intentDV_ib[sample(1:nrow(ctrl.df),i,replace=TRUE)],na.rm=TRUE), median(trt.df$intentDV_ib[sample(1:nrow(trt.df),i,replace=TRUE)], na.rm=TRUE))
		}
		ig.ib[h,i] <- mean(temp.df$DV[which(temp.df$trt==1)],na.rm=TRUE) - mean(temp.df$DV[which(temp.df$trt==0)],na.rm=TRUE)
	}
}
	
ig.ib2 <- as.data.frame(ig.ib)
colnames(ig.ib2) <- c(1:20)
ig.ib2$id <- 1:nrow(ig.ib2)

ig.ib.df <- pivot_longer(ig.ib2, cols=1:20)
ig.ib.df$name <- as.numeric(ig.ib.df$name)


#Reactive devaluation group size analysis for individual (generate 1000 groups with replacement (all group members must be in the same treatment condition), each of size N (where N=1-20); calcualte ATE, and resample 1000 times; repeat for each group size)
set.seed(43215)
B <- 1500
ig.rd <- matrix(NA, nrow=B, ncol=20)
ctrl.df <- dat.i[which(dat.i$chinaTrt_rd==0),] #Subset of observations in control condition
trt.df <- dat.i[which(dat.i$chinaTrt_rd==1),] #Subset of observations in treatment condition
for (h in 1:B){
	for (i in 1:20){
		temp.df <- data.frame(trt=rbinom(1000,1,0.5), DV=NA)	
		for (j in 1:length(temp.df$trt)){
			temp.df$DV[j] <- ifelse(temp.df$trt[j]==0,median(ctrl.df$supportDV_rd[sample(1:nrow(ctrl.df),i,replace=TRUE)],na.rm=TRUE), median(trt.df$supportDV_rd[sample(1:nrow(trt.df),i,replace=TRUE)], na.rm=TRUE))
		}
		ig.rd[h,i] <- mean(temp.df$DV[which(temp.df$trt==1)],na.rm=TRUE) - mean(temp.df$DV[which(temp.df$trt==0)],na.rm=TRUE)
	}
}
	
ig.rd2 <- as.data.frame(ig.rd)
colnames(ig.rd2) <- c(1:20)
ig.rd2$id <- 1:nrow(ig.rd2)

ig.rd.df <- pivot_longer(ig.rd2, cols=1:20)
ig.rd.df$name <- as.numeric(ig.rd.df$name)

#Merge together
df.all <- rbind(ig.df, ig.ib.df, ig.rd.df)
df.all$module <- as.factor(rep(c("Prospect theory", "Intentionality bias", "Reactive devaluation"),each=30000))
df.all$module <- relevel(df.all$module, ref="Prospect theory")

dev.new(height=5,width=10)
ggplot(df.all, aes(x=name, y=value, group=name)) + geom_boxplot() + xlab("Simulated group size") + ylab("ATE") + theme_bw() + facet_wrap(~module, scales="free") + geom_hline(aes(yintercept=0), lty=3)


##### Group participation analysis

#Descriptive statistics in text:

## What proportion of group members participated more than once in each deliberation?

#Horizontal
prop.table(table(dat.f$timesSpoken_pt>1)) #76% of group members participate more than once in the chat
prop.table(table(dat.f$timesSpoken_ib>1)) #73% of group members participate more than once in the chat
prop.table(table(dat.f$timesSpoken_rd>1)) #73% of group members participate more than once in the chat

#Hierarchical
#Proportion of leaders who speak at least once in a given module 
prop.table(table(dat.h$timesSpoken_pt[which(dat.h$role_pt=="Leader")] > 1)) #86% of group leaders participate more than once in the chat
prop.table(table(dat.h$timesSpoken_ib[which(dat.h$role_ib=="Leader")] > 1)) #82% of group leaders participate more than once in the chat
prop.table(table(dat.h$timesSpoken_rd[which(dat.h$role_rd=="Leader")] > 1)) #77% of group leaders participate more than once in the chat

# Proportion of advisers who speak at least once in a given module
prop.table(table(dat.h$timesSpoken_pt[which(dat.h$role_pt!="Leader")] > 1)) #80% of group advisers participate more than once in the chat
prop.table(table(dat.h$timesSpoken_ib[which(dat.h$role_ib!="Leader")] > 1)) #76% of group advisers participate more than once in the chat
prop.table(table(dat.h$timesSpoken_rd[which(dat.h$role_rd!="Leader")] > 1)) #74% of group advisers participate more than once in the chat


### Table 4.9: Little evidence of heterogeneous effects by participation


pt.f <- lm(medDV_pt ~ lossTrt_pt*timesSpokenPC_pt_g, data=dat.f.1.pt)
ib.f <- lm(medDV_ib ~ fatalTrt_ib*timesSpokenPC_ib_g, data=dat.f.1.ib)
rd.f <- lm(medDV_rd ~ chinaTrt_rd*timesSpokenPC_rd_g, data=dat.f.1.rd)

pt.h <- lm(riskyDV_pt ~ lossTrt_pt*timesSpokenPC_pt_g, data=dat.h.1.pt)
ib.h <- lm(intentDV_ib ~ fatalTrt_ib*timesSpokenPC_ib_g, data=dat.h.1.ib)
rd.h <- lm(supportDV_rd ~ chinaTrt_rd*timesSpokenPC_rd_g, data=dat.h.1.rd)

pt.i <- lm(riskyDV_pt ~ lossTrt_pt*wordsSpoken_pt, data=dat.i)
ib.i <- lm(intentDV_ib ~ fatalTrt_ib*wordsSpoken_ib, data=dat.i)
rd.i <- lm(supportDV_rd ~ chinaTrt_rd*wordsSpoken_rd, data=dat.i)

#To make output table relatively painless, rename coefficients
names(pt.f$coefficients) <- c("(Intercept)", "Treatment", "Group participation", "Trt x Group participation")
names(pt.h$coefficients) <- c("(Intercept)", "Treatment", "Group participation", "Trt x Group participation")
names(ib.f$coefficients) <- c("(Intercept)", "Treatment", "Group participation", "Trt x Group participation")
names(ib.h$coefficients) <- c("(Intercept)", "Treatment", "Group participation", "Trt x Group participation")
names(rd.f$coefficients) <- c("(Intercept)", "Treatment", "Group participation", "Trt x Group participation")
names(rd.h$coefficients) <- c("(Intercept)", "Treatment", "Group participation", "Trt x Group participation")
names(pt.i$coefficients) <- c("(Intercept)", "Treatment", "Individual participation", "Trt x Indiv participation")
names(ib.i$coefficients) <- c("(Intercept)", "Treatment", "Individual participation", "Trt x Indiv participation")
names(rd.i$coefficients) <- c("(Intercept)", "Treatment", "Individual participation", "Trt x Indiv participation")

stargazer(pt.f, pt.h, pt.i, ib.f, ib.h, ib.i, rd.f, rd.h, rd.i, title="Little evidence of heterogeneous effects by participation", style="apsr", digits=3, dep.var.labels.include=FALSE, omit.stat=c("LL", "ser", "f"), column.labels=c("Prospect theory", "Intentionality bias", "Reactive devaluation"), column.separate=c(3,3,3), add.lines=list(c("Group condition","Horizontal", "Hierarchical", "Individual","Horizontal", "Hierarchical", "Individual","Horizontal", "Hierarchical", "Individual")))